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Abstract: 

Formal Laurent-Puiseux series are important in many branches of mathemat- 
ics. This paper presents a Mathematica implementation of algorithms devel- 
oped by the author for converting between certain classes of functions and 
their equivalent representing series. The package PowerSeries handles func- 
tions of rational, exponential, and hypergeometric type, and enables the user 
to reproduce most of the results of Hansen's extensive table of series. Subal- 
gorithms of independent significance generate differential equations satisfied 
by a given function and recurrence equations satisfied by a given sequence. 



Scope of the Algorithms 

A common problem in mathematics is to convert an expression involving el- 
ementary or special functions into its corresponding formal Laurent-Puiseux 
series of the form 

oo 

f(x) = ^x k/n . (1) 

k=ko 

Expressions created from algebraic operations on series, such as addition, 
multiplication, division, and substitution, can be handled by finite algo- 
rithms if one truncates the resulting series. These algorithms are imple- 
mented in Mathematical Series command. For example: 

In[l]:= Series [Sin [x] Exp [x] , {x, 0, 5}] 

3 5 
2 x x 6 

0ut[l]= x + x + — + 0[x] 

3 30 

It is usually much more difficult to find the exact formal result, that is, an 
explicit formula for the coefficients a^. 
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This paper presents the package PowerSeries, a Mathematica imple- 
mentation of algorithms developed by the author (|J, j|, ||) for comput- 
ing Laurent-Puiseux series for three types of functions: functions of rational 
type, which are rational or have a rational derivative of some order; functions 
of exponential type, which satisfy a homogeneous linear differential equa- 
tion with constant coefficients; and functions of hypergeometric type, which 
have a representation (|l]) with coefficients satisfying a recurrence equation 
of the form 

ak+m = R{k) a k for k > k (2) 
a,k = for k = ko, ko + 1, • • • , &o + m ~ 1 

for some rational function R, A k £ C, Ak ^ 0, and some positive integer m, 
called the symmetry number of (the given representation of) the function. 

To find the Laurent-Puiseux expansion of a function f(x) about the 
point xq, one uses the function call PowerSeries [/(x) , {x , Xo}] • The 
package also extends the built-in Series command in the case that its second 
argument is a list of length two: 

In [2]:= << PowerSeries .m 

ps-info: PowerSeries version 1.02, Mar 07, 1994 

In [3]:= Series [Sin [x] Exp [x] , {x, 0}] 

k/2 k k Pi 

2 x Sin[ ] 

4 

Out [3]= Sum[ , {k, 0, Infinity}] 

k! 

Our algorithms also handle the converse problem of finding a "closed 
form" representation of a given Laurent-Puiseux series (that is, the gen- 
erating function of the sequence a^). The standard Mathematica pack- 
ages Algebra' SymbolicSum' and DiscreteMath'RSolve ' also deal with 
this problem, using different approaches. We will discuss the connection 
of our package with these packages and other Mathematica functions at the 
end of this introduction. 

The most interesting case is formed by the functions of hypergeomet- 
ric type, which include almost all transcendental elementary functions like 
x"n, Exp, Log, and the trigonometric and inverse trigonometric functions; all 
kinds of special functions like the Airy functions, the Bessel functions, the in- 
tegral functions Erf, ExpIntegralEi, Coslntegral, and Sinlntegral, the 
associated Legendre functions LegendreP and LegendreQ, the hypergeomet- 
ric functions HypergeometricU, HypergeometricOFl, HypergeometriclFl, 
as well as Hypergeometric2Fl, the orthogonal polynomials JacobiP, 
GegenbauerC, ChebyshevT, ChebyshevU, LegendreP, LaguerreL, and 
Hermit eH; and many more functions. 
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Our algorithms also handle the following special functions. Since they 
are not built in to Mathematica, we have included the necessary defini- 
tions in our package: the Bateman functions Bateman[n,x] ([Q], (13.6), and 
H); the Hankel functions Hankell[n,x] and Hankel2[n,x] (JTJ, (9.1)); 
the Kummer and Whittaker functions KummerM[a,b,x] , KummerU[a,b,x] , 
WhittakerM[a,b,x] , WhittakerW[a,b,x] (0), (13.1)); the Struve func- 
tions StruveH [n,x] and StruveL [n,x] (Q, Chapter 12); the hypergeo- 
metric functions HypergeometriclFO [a,x] , Hypergeometric2F0 [a,b,x] ; 
the iterated integrals of the complementary error function Erf c [n , x] (|M , 
(7.2)); the Abramowitz functions Abramowitz [n,x] ([ffl], (27.5)); and the 
repeated integrals of the normal probability integral Normallntegral [n,x] 
(|], (26.2.41)). Informat ions about these functions can be obtained by using 
Mathematical ? feature after the package is loaded. 

Our algorithms depend strongly on the fact that any function of hyper- 
geometric type satisfies a simple differential equation, i. e., a homogeneous 
linear differential equation with polynomial coefficients (see [||, Theorems 
2.1, and 8.1). One of our main algorithms finds this differential equation for 
(practically) any expression for which such an equation exists. 

On the other hand, it is similarly important that almost every function 
that we may write down satisfies a simple differential equation. All functions 
that can be constructed algebraically from the functions mentioned above 
by addition, multiplication, and composition with rational functions and 
rational powers, satisfy such a differential equation ( [fTCf] , Q). 

The conversion of functions to and from their representing Laurent- 
Puiseux expansions is described in detail in (||, Q). Here, we give a brief 
description and some examples for each of the types of functions considered. 

To find the Laurent-Puiseux expansion of an expression f(x), the func- 
tion PowerSeries performs the following steps. First, a homogeneous linear 
differential equation with polynomial coefficients for / is generated. (If no 
low-order differential equation exists, this may take some time.) This differ- 
ential equation is converted to an equivalent recurrence equation for at, an 
easy and fast procedure. If the recurrence equation is of the hypergeomet- 
ric type, the coefficients can be calculated explicitly from a finite number 
of initial values. (Finding the initial values may require the calculation of 
limits. This procedure may be slow, or may fail.) 

As an example of a function of rational type, let's take f(x) = arctanx. 
Here, f'(x) = j^-p is rational, so we apply a special algorithm for the ratio- 
nal case, and integrate. First, the complex partial fraction decomposition is 
derived: 

1 1 / 1 1 \ 

/ (X) ~ TT^2 - 2 \TTix' + ~Y^i~x) ' 

oo 

from which the coefficients of the derivative F'(x) = bkX can be 

fc=o 

deduced: 
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Finally, after an integration, we get 

°° i k (l + (-l) k , 

To illustrate the hypergeometric case, we again take f{x) = arctanx. Wc 
have f'(x) = and f"(x) = — jp^ji i so we arrive at the differential 

equation 

(l + x 2 )f"(x) + 2xf'(x) = 0. 
This differential equation can be converted to the recurrence equation 

(k + 1) (k + 2) a k+2 + k (k + 1) a k = (3) 

by the rule-based function detore [DE,F,x,a,k] , which takes as input the 
left hand side of the differential equation for the function F [x] and computes 
the left hand side of the recurrence equation for the coefficients a [k] : 

detore [g_ + h_, F_, x_, a_, k_] := 

detore [g, F, x, a, k] + detore [h, F, x, a, k] 
detore[c_ * g_, F_, x_, a_, k_]:= c * detorelg, F, x, a, k] /; 

FreeQlc, x] && FreeQ[c, F] 
detore [Derivative [k0_] [F_] [x_] , F_, x_, a_, k_]:= 

Pochhammer [k+1, kO] * a[k+k0] 
detore [F_ [x_] , F_, x_, a_, k_]:= a[k] 

detore [x_~j_ . * Derivative [k0_] [F_] [x_] , F_, x_, a_, k_]:= 

Pochhammer [k+l-j , kO] * a[k+k0-j] 
detore [x_~j_ . * F_ [x_] , F_, x_, a_, k_]:= a[k-j] 

In [4]:= Simplify [detore [2*x*f ' [x] + f " [x] + x~2*f"[x], f, x, a, k] ] 

Out [4]= (1 + k) (k a[k] + 2 a [2 + k] + k a [2 + k] ) 

As a last step, this recurrence equation is solved, resulting in the represen- 
tation 



/(*) = E 



Note that in most cases the hypergeometric-type result, if applicable, is the 
simplest one. 

Finally, let us take f(x) = sinxe^, the function whose series was calcu- 
lated above. Then f'(x) = e^suix + cosx^J and f"(x) = 2e :E cosa;. If we 
look for an equation of the form 

/" + A 1 f' + Aof = 2e x cos x + A x e x (sin x + cos xj + A e x sin x = , 

we may expand and recombine rationally dependent expressions to get 

/" + A 1 f + Aof = (2 + A i y x cos x+ (A 1 +A o y x sin x = 0. 
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This equation must be valid for all x, so the coefficients of e x cos x and 
e x sin x must vanish. Hence A\ = —2, Aq = 2, and we have the differential 
equation for /: 

/" - 2/' + 2/ = . 

This differential equation has constant coefficients, so / is of exponential 
type. For b n = n\ a n , we have the recurrence equation 



b n+2 - 2b n+ i + 2b n = 

with initial values 

bo = ao = /(0) = and b\ = a\ 

which generates the solution 

b n l(l + i) n -(I-*)' 



f'(0) 
1! 



nl n\ 2i 
—Jm (1 + i) n = irlm (V2e^) n 



nl nl 

1 n mr 
— 22 sm — 
n\ 4 



Hence 



f (rc) = > —22 sm — x 
w ^ n 4 



n 



n=0 

We will present more examples to demonstrate our Mathematica implemen- 
tation later in the paper. Additional examples are also given in Q. 

To find the generating function of a given sequence, the process can be 
reversed. The function Convert provides the following procedure to cal- 

oo 

culate the closed- form representation of • First, a homogeneous 

k=ko 

linear recurrence equation with polynomial coefficients for at is generated. 
(If no low order recurrence equation exists, this may take some time.) This 
recurrence equation is converted to an equivalent differential equation for 
the generating function, also an easy and fast procedure. Finally, the differ- 
ential equation is solved using appropriate initial values. (This step depends 
heavily on Mathematica' s built-in differential equation solver DSolve. It of- 
ten takes a lot of time to solve the differential equation, or may also fail. 
Moreover, the calculation of limits is needed to find the initial values. This 
part of the procedure may be slow or fail, as well.) 

As an example, we try to find the generating function 

of the sequence = 2fc+1 . Taking a shift, we get the relation ak+i = 2 fc+3 
so that the recurrence equation 

(2k + 3)a k+1 + (2k + l)a k = 
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is established for k > 0. If we multiply this equation by the factor {k + 1), 
the resulting recurrence equation 

(k + l)(2k + 3K+i + (k + l)(2k + l)a k = (4) 

holds for all integers k. The back-substitution can be done by the Mathe- 
matica procedure retode [re ,n,f ,x] , defined by these rules: 

theta[f_, x_] := x D [f , x] 

retode [eql_ + eq2_, k_ , f_, x_]:= 

retode[eql, k, f, x] + retode[eq2, k, f, x] 
retode[c_ * eq_, k_, f_, x_]:= c * retode[eq, k, f, x] /; 

(FreeQ[c, k] && FreeQ[c, f] kk FreeQ[c, x] ) 
retode [a [k_ + m_.], k_ , f_, x_]:= f[x]/x~m 
retode[k~j_. * a[k_ + m_.], k_, f_, x_]:= 

theta[retode [k~ (j-l) * a[k+m] , k, f, x] , x] 
retode [p_ * a[k_ + m_.], k_, f_, x_]:= 

retode [Expand [p * a[k+m]], k, f, x] /; PolynomialQ [p] 

Applied to the left hand side of (0), retode gives: 

In[5]:= Simplify [retode [ (k+1) (2k+3) a [k+1] + (k+1) (2k+l) a [k] , k, f, x] ] 
Out [5] = f [x] + 3 f ' [x] + 5 x f ' [x] + 2 x f " [x] + 
2 

> 2 x f " [x] 

The initial value problem 

(2x + 2x 2 )f" + (3 + 5x)f + / = , 

/(0) = a = l, f(0)=ai = -i, 

has the solution 

arctan y/x 



Note that Mathematica Version 2.0 gives the correct result: 

In [6]:= Convert[Sum[(-l)~k/(2k+l)*x~k, {k, 0, Infinity}], x] 

ArcTan [Sqrt [x] ] 

Out [6]= 

Sqrt [x] 

whereas the DSolve command of Mathematica Version 2.2 gives a wrong 
result. 

Now we want to give some remarks on the connection of our package 
with several standard Mathematica packages. 
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In several instances, PowerSeries depends on the calculation of limits, so 
it may be helpful to load the package Calculus ' Limit ' . Since the Convert 
function of PowerSeries depends heavily on the DSolve command, it may 
also be helpful to load Calculus 'DSolve ' (Version 2.2). In this case, the 
incorrect result mentioned above does not occur. 

The aim of the package Algebra' SymbolicSum' is to give closed form 
representations of sums, in particular, sums of the generating function type. 
Thus, it is not surprising that many of the conversions that can be obtained 
by our Convert function are also obtained by the SymbolicSum package. 
Among the expressions that cannot be handled by the SymbolicSum package 
are the sums conv[4] and conv[5] defined below, and, in general, sums 
involving special functions. Typically, our implementation is slower as it 
always uses DSolve, which takes a lot of time. The advantage is that the 
procedure forms an algorithmic treatment of the problem, and is not based 
on heuristics. 

The purpose of the package DiscreteMath'RSolve ' is to solve recur- 
rence equations. It uses generating function techniques and defines the 
function GeneratingFunction to find the generating function of a recur- 
rence equation. The procedure uses heuristic techniques, and in particular 
treats hypergeometric recurrence equations (equation (||) with m = 1), but 
the general case (to > 1) is not implemented. Similar restrictions apply to 
the RSolve command. 

We give an example: Our function SimpleRE generates the recurrence 
equation (Q) for the coefficients of the series representation of a function: 

In [7]:= re=SimpleRE [ArcTan [x] , x] 

Out [7]= k (1 + k) a[k] + (1 + k) (2 + k) a[2 + k] ==0 

We load the package and try to solve the recurrence equation: 

In [8]:= « DiscreteMath'RSolve' 

In [9] := RSolve [{re, a[0]==0, a[l]==l>, a[k] , k] 

DSolve : : dnim: 

Built-in procedures cannot solve this differential 
equation. 

Out [9]= RSolve [{re, a[0] == 0, a[l] == 1}, a[k], k] 

(This is the Mathematica Version 2.2 result. In Version 2.0, another error 
message occurs.) We see that RSolve tries to solve the differential equa- 
tion for the generating function of the sequence a k given by the recurrence 
equation 

k(l + k)a k + (l + k)(2 + k)a k+2 = 0. (5) 



7 



This problem is too hard. In our treatment, the recurrence equation is 
replaced by the simpler equation 

k (1 + k) a k + (1 + k) (2 + k) a k+1 = 

whose generating function is easier to find. Further, our function PowerSeries 
solves the recurrence equation (||) by a direct application of the hypergeo- 
metric theory, and splits the series into two parts (odd and even), one of 
which disappears as a result of a vanishing initial value: 

In [10]:= Series [ArcTan[x] , {x, 0}] 

k 1 + 2 k 
(-1) x 

Out [10]= Sum[ , {k, 0, Infinity}] 

1 + 2 k 



Generation of Differential Equations 

In this section, we present some examples of the function SimpleDE, which 
generates the homogeneous linear differential equation of least order satisfied 
by a given function. 

In [11]:= de[l] = SimpleDE [x~n * E~ (alpha x) , x] 

Out [11]= (-n - alpha x) F [x] + x F' [x] == 

In[12]:= de[2] = SimpleDE [((l+x)/(l-x))~n, x] 

Out [12]= 2 n F[x] + (-1 + x) (1 + x) F' [x] == 

In [13]:= de[3] = SimpleDE [ArcSin [x"5] , x] 

10 11 
Out [13]= (4 + x ) F'[x] + (-x + x ) F" [x] == 

In [14]:= de[4] = SimpleDE [ArcSin [x] , x] 

2 

Out [14]= x F'[x] + (-1 + x ) F"[x] == 

As noted above, we can create new functions satisfying this kind of differ- 
ential equation by addition, multiplication, and by the composition with 
rational functions and rational powers: 

In [15]:= de[5] = SimpleDE [ArcSin [x] "3 , x] 

2 
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Out [15]= x F' [x] + (-4 + 7 x ) F" [x] + 

2 (3) 2 2 (4) 

> 6 x (-1 + x ) F [x] + (-1 + x ) F [x] == 

In [16]:= de[6] = SimpleDE [Sin [x] "5, x, 6] 

(4) 

Out [16]= 225 F[x] + 259 F"[x] + 35 F [x] + 
(6) 

> F [x] == 

In [17]:= de[7] = SimpleDE [Exp [alpha x] * Sin [beta x] , x] 
2 2 

Out [17]= (alpha + beta ) F[x] - 2 alpha F' [x] + 

> F" [x] == 

By default, SimpleDE searches for a simple differential equation of order up 
to five. This setting can be changed by entering the order as a third argu- 
ment. In case one knows a priori that a differential equation of the given type 
exists, one may choose the order Infinity. Note that there is no simple dif- 
ferential equation of order less than five for sin 5 x, so SimpleDE [Sin [x] "5 , x] 
fails. 

Our implementation handles orthogonal polynomials and special func- 
tions using the general method described in Q. We can derive the usual 
differential equations: 

In [18]:= de[8] = SimpleDE [LaguerreL [n, x] , x] 

Out [18]= n F[x] + (1 - x) F'[x] + x F" [x] ==0 

In [19]:= de[9] = SimpleDE [ChebyshevT[n, x] , x] 

2 2 
Out [19]= -(n F[x] ) + x F' [x] + (-1 + x ) F" [x] == 

In [20]:= de[10] = SimpleDE [BesselY [n, x] , x] 

2 2 2 
Out [20]= (-n + x ) F[x] + x F' [x] + x F" [x] == 

In [21]:= de[ll] = SimpleDE [AiryAi [x] , x] 

Out [21]= -(x F[x]) + F"[x] == 

By the statement above, we know that much more complicated functions 
satisfy differential equations: 
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In [22]:= de[12] = S imp leDE [Exp [alpha x] * Bessell[n, x] , x] 

2 2 2 2 

Out [22]= (-n - alpha x - x + alpha x ) F[x] + 

2 2 

> (x - 2 alpha x ) F' [x] + x F" [x] == 

In[23]:= de[13] = SimpleDE [Sin [m x] * BesselJ[n, x] , x] 

2 4 62 22 22 
Out [23]= (n -5n +4n -x -2m x +22n x - 

222 42 242 4 

> 10mnx-12nx+12mnx-x + 

24 44 24 224 

> 6mx-5mx+12nx-24mnx + 

424 6 26 46 

> 12mnx-4x+12mx-12mx + 

6 6 

> 4 m x ) F[x] + 

3 23 23 223 

> (-8 x -4m x +8n x +40m n x - 

5 4 5 

> 8x+8mx)F'[x] + 

2 22 42 4 24 

> (-2 x + 10 n x -8n x +2x -6m x + 

2 4 6 4 6 

> 16 n x - 8x + 8m x ) F" [x] + 

3 2 3 5 2 5 (3) 

> (-4 x + 16 n x - 8 x + 8 m x ) F [x] + 

4 2 2 2 2 (4) 

> x (-1 + 4 n - 4 x + 4 m x ) F [x] == 

In[24]:= de[14] = SimpleDE [LegendreP[n, x]~2, x] 

2 

Out [24]= (-4 n x - 4 n x) F[x] + 

2 2 2 2 2 

> (-2 +4n+4n +6x -4nx -4n x) 
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2 

> F' [x] + 6 x (-1 + x ) F' ' [x] + 



2 2 (3) 

> (-1 + x ) F [x] == 

(The calculation of these last two results is very time consuming!) 
The functions 

where iffi denotes the generalized Laguerre polynomial, were studied by the 
author in Q. For these functions, we immediately deduce the differential 
equation: 

In[25]:= de[15] = SimpleDE [E~ (-x) * LaguerreL[n, alpha, 2x] , x] 
Out [25]= (1 + alpha + 2 n - x) F [x] + 

> (1 + alpha) F' [x] + x F" [x] == 

We may observe that the coefficient of F' in the differential equation vanishes 
only in the case a = — 1, which is studied in M. 

Conversion of Simple Differential Equations to Re- 
currence Equations 

The function DEtoRE converts a simple differential equation into a recurrence 
equation for the coefficients of a corresponding Laurent-Puiseux expansion. 
The output is given in terms of the symbols a and k representing the coeffi- 
cient sequence at-- As examples, we convert some of the differential equations 
that were derived above. 
The example 

In [26]:= re[l] = DEtoRE [de [2] , F, x] 

Out [26]= k a[k] + 2 n a[l + k] - (2 + k) a [2 + k] ==0 

shows that the function (x^f ) is not of hyper geometric type, whereas the 
example 

In [27]:= re [2] = DEtoRE [de [3] , F, x] 
2 

Out [27]= k a[k] - (5 + k) (10 + k) a [10 + k] ==0 

shows that the function arcsin(x 5 ) is of hypergeometric type with symmetry 
number m = 10. Similarly, the Airy function is of hypergeometric type with 
symmetry number m = 3: 
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In [28]:= re [3] = DEtoRE [de [11] , F, x] 

Out [28]= -a[k] + (2 + k) (3 + k) a [3 + k] ==0 

Among the functions of the form e ax I n (x), where a £ H and / is the 
modified Bessel function of the first kind, only the cases a = and a = ±1 
represent functions of hyper geometric type, as is seen from the recurrence 
equation 

In [29]:= re [4] = DEtoRE [de [12] , F, x] 
Out [29]= (-1 + alpha) (1 + alpha) a[k] - 

> alpha (3 + 2 k) a[l + k] + 

> (2 + k - n) (2 + k + n) a [2 + k] == 

Similarly, the functions of the form sin (mx) J m (x), where (m £ II) and J 
is the Bessel function of the first kind, are of hyper geometric type only for 
the cases m = and m = ±1, as is seen from the recurrence equation 

In[30]:= re [5] = DEtoRE [de [13] , F, x] 

3 3 
Out [30]= 4 (-1 + m) (1 + m) a[k] + 

> (-1 + m) (1 + m) 

2 2 2 2 2 

> (33 + 32 k + 8 k +27m +32 km + 8k m - 

2 2 2 

> 12 n + 12 m n ) a [2 + k] + 

2 3 4 2 

> (-297 - 402 k - 210 k - 48 k -4k + 198 m + 

2 2 2 3 2 4 2 

> 362 k m + 206 km+48km+4km+ 

2 2 2 2 2 2 

> 246 n + 120 k n + 16 k n + 150 m n + 

2 2 4 2 4 

> 40 k m n - 12 n + 12 m n ) a [4 + k] + 

> (5 + k - n) (6 + k - n) (5 + k + n) (6 + k + n) 

> (-1 + 2 n) (1 + 2 n) a [6 + k] ==0 
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The function SimpleRE combines the two steps SimpleDE and DEtoRE to 
calculate the recurrence equation of a Laurent-Puiseux series representation 
of a function. 



Calculation of Laurent-Puiseux Series 

In this section, we present some of the results of the procedure PowerSeries. 
Obviously, the program can handle the standard power series representations 
of the elementary functions: 

In [31]:= ps[l] = PowerSeries [E~x, x] 

k 

x 

Out [31]= Sum[— , {k, 0, Infinity}] 
k! 



In [32]:= ps [2] = PowerSeries [Sin [x] , x] 

k 1 + 2 k 
(-1) x 

Out [32]= Sum[ , {k, 0, Infinity}] 

(1 + 2 k) ! 



In [33]:= ps [3] = PowerSeries [ArcSin [x] , x] 

1 k 1 + 2 k 2 
(-) x (2k)! 
4 

Out [33]= Sum[ , {k, 0, Infinity}] 

2 

k! (1 + 2 k) ! 



Here is a series expanded about the point x = 1: 
In [34]:= ps [4] = PowerSeries [Log [x] , {x, 1}] 

k 1 + k 

(-1) (-1 + x) 

Out [34]= Sum[ , {k, 0, Infinity}] 

1 + k 

It may come as a surprise that the following functions are of hypergeometric 
type: 

In [35] := ps [5] = PowerSeries [Exp [ArcSin [x] ] , x] 
k 2 k 5 2 
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4 x Product [- - 2 jj + jj , -[jj, k>] 
4 

Out [35]= Sum[ , 

(2 k)! 

> {k, 0, Infinity}] + 

k 1 + 2 k 1 2 

4 x Product [- - jj + jj , {jj, k}] 

2 

> Sum [ , 

(1 + 2 k) ! 

> {k, 0, Infinity}] 

In [36] := ps [6] = PowerSeries [Exp [ArcSinh [x] ] , x] 

1 k 2 k 
(-(-)) x (2k)! 
4 

Out [36]= x + Sum[ , {k, 0, Infinity}] 

2 

(1 - 2 k) k! 

If one is interested in the intermediate calculations, one may give the input psprint, 
which sets the global variable PSPrintMes sages to True: 

In [37] := psprint 

In [38]:= ps [7] = PowerSeries [E~x - 2 E~(-x/2) Cos [Sqrt [3] x/2 - Pi/3], x] 
ps-info: PowerSeries version 1.02, Mar 07, 1994 
ps-info: 3 step(s) for DE: 
(3) 

-F[x] + F [x] == 
ps-info: RE for all k >= : 







a [3 


+ k] 


ps- 


•inf o : 


function 


ps- 


info : 


a[0] 


= 


ps- 


info : 


a[l] 


= 
3 


ps- 


info : 


a [2] 


2 


ps- 


info : 


a [3] 


= 


ps- 


info : 


a [4] 


= 



2 + 3 k 

9 (1 + k) x 

Out [38]= Sum[ , {k, 0, Infinity}] 

(3 + 3 k) ! 

Note that the differential equation in this example is extremely simple. 
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We now suppress intermediate results using the command nopsprint, which 
sets PSPrintMessages to False. 

In [39] := nopsprint 

oo 

Here is a well-known closed formula for the generating function f(x) — o,kX k 

fc=0 

of the Fibonacci numbers a n defined by the recurrence 

a n+ \ =«„ + a n _ i , ao = , oi = 1 . 

In [40]:= ps[8] = PowerSeries [x/(l - x - x~2) , x] 

k 1 k -2 k k 

(2 ( ) _ ( ) ) x 

-1 + Sqrt[5] 1 + Sqrt[5] 

Out [40]= Sum[ , 

Sqrt [5] 

> {k, 0, Infinity}] 

Moreover, we can use SimpleRE to generate the Fibonacci recursion: 
In [41]:= re [6] = SimpleRE [x/(l - x - x~2) , x, a, k] 

Out [41]= (1 + k) a[k] + (1 + k) a [1 + k] - 

> (1 + k) a[2 + k] ==0 

Note that the common factor (1 + fc) ensures that the recurrence equation holds for 
all integers k. 

Our implementation covers functions that correspond to hypergeometric type 
Laurent-Puiseux series, rather than just power series. For examples, see Q]. 

Here, we want to return to the case of special functions. In the previous section, 
we discovered that the functions e x I n {x) and sin (x) J n (x) are of hypergeometric 
type. For n = and 1, we get the following series representations: 

In [42]:= ps [9] = PowerSeries [E~x * Bessell[0, x] , x] 

1 k k 
(-) x (2k)! 
2 

Out [42]= Sum[ , {k, 0, Infinity}] 

3 

k! 

In [43]:= ps [10] = PowerSeries [E~x * Bessell[l, x] , x] 

lk 1 + k 
(-) x (1 + 2 k) ! 

2 

Out [43]= Sum[ , {k, 0, Infinity}] 

2 

k! (2 + k) ! 
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In [44]:= ps[ll] = PowerSeries [Sin [x] * BesselJ[0, x] , x] 

1 k 1 + 2 k 
(-(-)) x (1 + 4 k) ! 

4 

Out [44]= Sum[ , 

2 

(2k)! (1 + 2 k) ! 

> {k, 0, Infinity}] 

In [45]:= ps [12] = PowerSeries [Sin [x] * BesselJ[l, x] , x] 

1 k 2 + 2 k 
(-(-)) x (3 + 4 k) ! 

4 

Out [45]= Sum[ , 

2 

2 (1 + 2 k) ! (3 + 2 k) ! 

> {k, 0, Infinity}] 

Many more functions constructed from the given special functions are of hyperge- 
ometric type (see, for example, Hansen's extensive table of series [^)), and their 
Laurent-Puiseux expansions can therefore be found by PowerSeries. 

Calculation of Generating Functions 

The function Convert calculates generating functions of series: 
In[46]:= conv[l] = Convert [Sum [(2k) ! /k ! "2 x~k, {k, 0, Infinity}], x] 
1 

Out [46]= 

Sqrt[l - 4 x] 

In [47]:= conv[2] = Convert [Sum [k ! "2/ (2k) ! x~k, {k, 0, Infinity}], x] 

Sqrt [x] 

4 Sqrt[x] ArcSin[ ] 

4 2 

Out [47]= + 

4 - x 3/2 
(4 - x) 

Again, with psprint we get information about the intermediate steps: 
In [48] := psprint 

In[49]:= conv[3] = Convert [Sum [k !/ (2k) ! x~k, {k, 0, Infinity}], x] 
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ps-info: PowerSeries version 1.02, Mar 07, 1994 
ps-info: 1 step(s) for RE. 

a[-l + k] + 2*(1 - 2*k)*a[k] == 
ps-info: DE: 

F[x] + (-2 + x) F'[x] + (-4*x) F"[x] == 
ps-info: Trying to solve DE . . . 
ps-info: DSolve 't encountered 
ps-info: DSolve computes 

C[3] + E-(x/4)*x-(l/2)* 

> (C[l] + (C[2]*(-2 - 

> E~ (x/4) *Pi~ (1/2) *x~ ( 1/2) *Erf [x~ (1/2) /2] ) ) / 

> (E~(x/4)*x~(l/2))) 
ps-info: expression rearranged: 

E~(x/4)*x~(l/2)*C[l] + C[3] + 

> E~ (x/4) *x~ (1/2) *C [2] *Erf [x~ ( 1/2) /2] 
ps-info: Calculation of initial values... 
ps-info: C [3] == 1 

ps-info: C[l] == 

ps-info: C [2] /Pi~ (1/2) == 1/2 

x/4 Sqrt [x] 

E Sqrt [Pi] Sqrt [x] Erf [ ] 

2 

Out [49]= 1 + 

2 

Here arc two generating function examples involving special functions. (Note 
that Convert can only be successful if the solution is representable by an elemen- 
tary function; otherwise, Mathematica will fail to solve the differential equation 
produced). 

In [50] := conv[4] = Convert [Sum [ChebyshevT[k,x] z~k, {k, 0, Infinity}], z] 

ps-info: PowerSeries version 1.02, Mar 07, 1994 
ps-info: 2 step(s) for RE. 

a [-2 + k] - 2*x*a[-l + k] + a[k] == 
ps-info: DE: 

2 F[z] + 4*(-x + z) F' [z] + (1 - 2*x*z + z~2)\ 

> F" [z] == 

ps-info: Trying to solve DE . . . 
ps-info: DSolve computes 

(z*C[l] - C[2] + (z*C[2])/C[l])/ 

> (-1 + 2*x*z - z~2) 
ps-info: expression rearranged: 

(z*C[l] + C[2])/(-l + 2*x*z - z~2) 



17 



ps-info: Calculation of initial values... 
ps-info: -C [2] == 1 

ps-info: (-2*C[1] - 4*x*C[2])/2 == x 

1 x z 

Out [50]= -( ) + 

2 2 

-l+2xz-z -l+2xz-z 

In[51]:= conv[5] = Convert [Sum [LaguerreL [k, a, x] z~k, {k, 0, Infinity}], z] 

ps-info: PowerSeries version 1.02, Mar 07, 1994 
ps-info: 2 step(s) for RE. 

(-1 + a + k)*a[-2 + k] + 

> (1 - a - 2*k + x)*a[-l + k] + k*a[k] == 
ps-info: DE: 

(-1 - a + x + z + a*z) F[z] + ((-1 + z)~2)\ 

> F' [z] == 

ps-info: Trying to solve DE . . . 
ps-info: DSolve computes 

E~(x/(-l + z))*(-l + z)-(-l - a)*C[l] 
ps-info: expression rearranged: 

E~(x/(-l + z))*(l - z)'(-l - a)*C[l] 
ps-info: Calculation of initial values... 
ps-info: C[l]/E~x == 1 

(x z)/(-l + z) -1 - a 

Out [51]= E (1 - z) 



Calculation of Recurrence Equations 

As described above, the function Convert works by generating a recurrence equa- 
tion for the given coefficients, coverting it to a differential equation, and then solv- 
ing. The procedure for generating the recurrence equation is implemented in the 
function FindRecursion. We give some examples: 

In[52] := re [7] = FindRecursion[(l + (-l)~n)/n, n] 
Out [52]= (2 - n) a [-2 + n] + n a[n] == 
In[53] := re [8] = FindRecursion [n + (-l)~n, n] 
Out [53]= (1 - 2 n) a[-2 + n] - 2 a[-l + n] + 
> (-3+2 n) a[n] == 

In[54] := re [9] = FindRecursion [(n + (-l)~n)/n~2, n] 
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2 3 

Out [54]= (4 - 12 n + 9 n - 2 n ) a [-2 + n] + 

2 2 3 

> (-2 + 4 n - 2 n ) a[-l + n] + (-3 n + 2 n ) a[n] 

> ==0 

In[55]:= re[10] = FindRecursion[l/(2n+l) ! , n] 

Out [55]= a[-l + n] - 2 n (1 + 2 n) a[n] == 

In [56] := ps [8] [[1]] 

k 1 k -2 k k 

(2 ( ) - ( ) ) x 

-1 + Sqrt[5] 1 + Sqrt[5] 

Out [56]= 

Sqrt [5] 



In[57]:= re[ll] = FindRecursion[7„, k] 
2 

0ut[57]= -(x a[-2 + k] ) - x a[-l + k] + a[k] == 

In[58]:= re[12] = FindRecursion [E~ (-x) LaguerreL [n, alpha, 2x] , n] 

Out [58]= (-1 + alpha + n) a [-2 + n] + 

> (1 - alpha - 2 n + 2 x) a[-l + n] + n a[n] == 
In[59] := re[13] = FindRecursion [n*LaguerreL [n, alpha, 2x] , n] 

2 

Out [59]= (1 - alpha - 2 n + alpha n + n ) a [-2 + n] + 

2 

> (-2 + 2 alpha + 5 n - alpha n-2n - 4 x + 

2 

> 2 n x) a[-l + n] + (2 - 3 n + n ) a[n] == 

In[60] := re[14] = FindRecursion [LaguerreL [n, alpha, 2x]/n, n] 

2 

Out [60]= (2-2 alpha - 3 n + alpha n + n ) 

> a [-2 + n] + (-1 + alpha + 3 n - alpha n - 

2 2 
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> 2 n - 2 x + 2 n x) a[-l + n] + n a[n] == 

In[61]:= re[15] = FindRecursion[LaguerreL[n, x]~2, n] 

2 3 2 

Out [61]= (4 - 12 i + 9 n - 2 n +4x-4nx+n x) 

2 3 

> a [-3 + n] + (-6 + 22 n - 21 n + 6 n - 14 x + 

2 2 2 3 

> 26 n x - 11 n x-7x +6nx - x ) 

2 3 

> a [-2 + n] + (2 - 10 n + 15 n - 6 n + 6 x - 

2 2 2 3 

> 18 e x + 11 n x+5x -6nx +x) 

2 3 2 

> a[-l + n] + (-3 n + 2 n - n x) a[n] == 
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Electronic Supplement 

The electronic supplement contains the package PowerSeries .m. 
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